# load packages
library(expss)
library(knitr)
library(sandwich)
library(lmtest)
library(plyr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(arm)
library(foreign)
library(lattice)
library(MASS)
library(stats)
library(base)
library(car)
library(xtable)
library(broom)
library(stringr)
library(grid)
library(gridExtra)
library(lme4)
library(stringr)
library(plm)
library(estimatr)
library(lmtest)
library(multiwayvcov)
library(modelr)
library(readr)
library(haven)
library(stargazer)
library(brms)
library(jtools)
library(sjPlot)
library(sjlabelled)
library(sjmisc)
library(ggeffects)
library(clubSandwich)
library(texreg)
library(fixest)
library(clusterSEs)
library(cregg)
library(cjoint)
library(knitr)
library(kableExtra)
library(magrittr)
library(logr)

rm(list = ls()) 

options(scipen=999)#turn off scientific notation or price of planes will be in scientific notation and the code for conjoint won't work.

multi <- read.csv("multi.csv") #load data


# group is the treatment variable. news is the reference level, then there are two other groups (generic info and specific conjoint info)
multi$group <- as.factor(multi$group)
multi$group <- relevel(multi$group, ref = "News")

multi$group <-factor(multi$group, levels = c("News",  "Specific Information", "Generic Information"))


# base models fairness is the DV 
mod1 <- lm(agree_1 ~ group, data=multi) 


summary(mod1)
pred1 <- ggpredict(mod1, terms = c("group")) 

#this is figure 1

pdf(file = "figure1.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(pred1) + 
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high, x=x, width = 0.1))+
  labs(x = "", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "DV: The company's decision to introduce the new technology is fair") + 
  theme(text = element_text(size = 18)) + ylim(1,5) + coord_flip() 
dev.off()

linearHypothesis(mod1, c("groupSpecific Information=groupGeneric Information"), white.adjust = "hc1")
#this tells us if generic and specific info coefficients are different from each other

## plot of policies (figure2)

summ <- multi %>%
  drop_na(group, policy_socialspend,policy_basicincome,policy_jobguarantee,
          policy_unskilled,policy_skilled,policy_trade,policy_reskill,policy_automationtax) %>%
  dplyr::group_by(group) %>%
  dplyr::summarize(socialsupport=mean(policy_socialspend, na.rm=T),
                   basicincome=mean(policy_basicincome, na.rm=T),
                   jobguarantee= mean(policy_jobguarantee, na.rm=T),
                   unskilled=mean(policy_unskilled, na.rm=T),
                   skilled=mean(policy_skilled, na.rm=T),
                   trade=mean(policy_trade, na.rm=T),
                   reskill=mean(policy_reskill, na.rm=T),
                   tax=mean(policy_automationtax, na.rm=T))

summ_long <- tidyr::gather(summ, policy_option, support, socialsupport: tax, 
                           factor_key=T)


pdf(file = "figure2.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
ggplot(summ_long, 
       aes(x=factor(policy_option), y=support, fill=group)) + 
  geom_bar(position="dodge", stat="identity") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 60, hjust=1, size=13),
        axis.text.y = element_text(size=14),
        axis.title.x = element_text(size=14),
        axis.title.y = element_text(size=14),
        plot.title = element_text(size=14, hjust=0.5, face="bold")) +
  ggtitle("") +
  ylim(0,5) +
  xlab("") +
  ylab("Mean Support") +
  scale_x_discrete(labels= c("Social support", "Basic Income", "Job Guarantee", "Unskilled", "Skilled",
                             "Trade", "Reskill", "Tax automation")) 
dev.off()

##############
# conjoint
###############
#select only variables we'll use
multi2 <- dplyr::filter(multi, group=="Specific Information") #this is the group that saw the conjoint
multi2 <- dplyr::select(multi2, numlsafter, numlsafter1, numlsafter2, numlsafter3, numhsafter, numhsafter1, 
                        numhsafter2, numhsafter3, wagelsafter, wagelsafter1, wagelsafter2, wagelsafter3, 
                        wagehsafter, wagehsafter1, wagehsafter2, wagehsafter3, priceafter, priceafter1, 
                        priceafter2, priceafter3, product, product1, product2, product3, scen1_agree_1,
                        scen2_agree_1, scen3_agree_1, scen4_agree_1,scen1_agree_2,
                        scen2_agree_2, scen3_agree_2, scen4_agree_2, ResponseId, UserLanguage, country,
                        scen1_know1, scen1_know2, female, college, LR, self_know_automation_cat, know_automation_correct,
                        screen1)

#rename variables 
multi2 <- dplyr::rename(multi2, numlsafter_1 = numlsafter, numlsafter_2 = numlsafter1, numlsafter_3 = numlsafter2,
                        numlsafter_4 = numlsafter3, numhsafter_1 = numhsafter, numhsafter_2 = numhsafter1, 
                        numhsafter_3 = numhsafter2, numhsafter_4 = numhsafter3, wagelsafter_1 = wagelsafter,
                        wagelsafter_2 = wagelsafter1, wagelsafter_3 = wagelsafter2, wagelsafter_4 = wagelsafter3, 
                        wagehsafter_1 = wagehsafter, wagehsafter_2 =wagehsafter1, wagehsafter_3 = wagehsafter2,
                        wagehsafter_4 = wagehsafter3, priceafter_1 = priceafter, priceafter_2 = priceafter1, 
                        priceafter_3 = priceafter2, priceafter_4 = priceafter3, product_1 = product, 
                        product_2 = product1, product_3 = product2, product_4 = product3, 
                        scenagree_1_1 = scen1_agree_1,
                        scenagree_2_1 = scen2_agree_1,  scenagree_3_1 = scen3_agree_1, 
                        scenagree_4_1 = scen4_agree_1, scenagree_1_2 = scen1_agree_2,
                        scenagree_2_2 = scen2_agree_2, scenagree_3_2 = scen3_agree_2, 
                        scenagree_4_2 = scen4_agree_2, scenknow_1_1 = scen1_know1 , scenknow_1_2 = scen1_know2)

#reshape dataset so we can analyze it
multi2_reshape <- multi2 %>% 
  #first step is gathering all the variables that we want to use as our attributes;
  #because we have named each task iteration in this format, using contains will get each conjoint-related variable
  gather(key = task, value = score, contains("_1"), contains("_2"), contains("_3"), contains("_4")) %>% 
  #we need to separate these variables out by characteristics; this can be done using `separate`, because variable were meaningfully named using underscores 
  separate(task, c("variable", "iteration", "index")) %>% 
  #creating an indicator of variables that were a part of an index
  dplyr::mutate(index = ifelse(is.na(index), "", index))  %>%
  #create a variable specific column
  unite("variable_index", c("variable", "index"), sep = "") %>% 
  #spreading out the variables so that each meaningful variable (or attribute) is a column
  tidyr::spread(key = variable_index, value = score) %>% 
  #lastly, because the data are a mix of upper and lower case, I use clean_names to make the variables names follow a uniform convention (this is optional, but recommended)
  janitor::clean_names() 

#create factor variables

multi2_reshape <- multi2_reshape %>%
  mutate(female = case_when(
    female== 0 ~"Male",
    female == 1 ~ "Female"
  ))

multi2_reshape <- multi2_reshape %>%
  mutate(college = case_when(
    college== 0 ~"No College",
    college == 1 ~ "College"
  ))

#makes sure all of our variables of interest are factors
multi2_reshape$female <- as.factor(multi2_reshape$female)
multi2_reshape$college <- as.factor(multi2_reshape$college)
multi2_reshape$know_automation_correct <- dplyr::recode(multi2_reshape$know_automation_correct,`0`="Low",`1`= "High")


multi2_reshape$scenagree1 <- as.numeric(as.character(multi2_reshape$scenagree1))
multi2_reshape$scenagree2 <- as.numeric(as.character(multi2_reshape$scenagree2))
multi2_reshape$numhsafter <- as.factor(as.character(multi2_reshape$numhsafter))
multi2_reshape$numlsafter <- as.factor(as.character(multi2_reshape$numlsafter))
multi2_reshape$priceafter <- as.factor(as.character(multi2_reshape$priceafter))
multi2_reshape$product <- as.factor(as.character(multi2_reshape$product))
multi2_reshape$wagehsafter <- as.factor(as.character(multi2_reshape$wagehsafter))
multi2_reshape$wagelsafter <- as.factor(as.character(multi2_reshape$wagelsafter))
multi2_reshape$self_know_automation_cat <- as.factor(as.character(multi2_reshape$self_know_automation_cat))

#recode attributes
multi2_reshape$priceafter <- ifelse(multi2_reshape$priceafter=="100000000" | 
                                      multi2_reshape$priceafter== "25" | 
                                      multi2_reshape$priceafter== "25000"| 
                                      multi2_reshape$priceafter== "600", "Same price", 
                                    ifelse(multi2_reshape$priceafter=="80000000" | 
                                             multi2_reshape$priceafter== "20" | 
                                             multi2_reshape$priceafter== "20000"| 
                                             multi2_reshape$priceafter== "480", "20% cheaper", 
                                           ifelse(multi2_reshape$priceafter=="50000000" | 
                                                    multi2_reshape$priceafter== "12.5" | 
                                                    multi2_reshape$priceafter== "12500"| 
                                                    multi2_reshape$priceafter== "300", "50% cheaper",  NA)))

multi2_reshape$product <- ifelse(multi2_reshape$product=="Vaccine" | 
                                   multi2_reshape$product=="Vaccin", "Vaccine",
                                 ifelse(multi2_reshape$product=="smartphone" | 
                                          multi2_reshape$product=="Téléphone intelligent", "Smartphone",
                                        ifelse(multi2_reshape$product=="Automobile" | 
                                                 multi2_reshape$product=="Car", "Car",
                                               ifelse(multi2_reshape$product=="Avion" | 
                                                        multi2_reshape$product=="Plane", "Plane",  NA))))

multi2_reshape$wagehsafter <- ifelse(multi2_reshape$wagehsafter=="£ 75,000" | 
                                       multi2_reshape$wagehsafter=="$125,000", "$125,000",
                                     ifelse(multi2_reshape$wagehsafter=="£ 90,000" | 
                                              multi2_reshape$wagehsafter=="$150,000", "$150,000",NA))

multi2_reshape$wagelsafter <- ifelse(multi2_reshape$wagelsafter=="£ 13,000" | 
                                       multi2_reshape$wagelsafter=="$20,000"| 
                                       multi2_reshape$wagelsafter=="$30,000",
                                     "$20,000",
                                     ifelse(multi2_reshape$wagelsafter=="£ 17,000" | 
                                              multi2_reshape$wagelsafter=="$25,000"| 
                                              multi2_reshape$wagelsafter=="$40,000",
                                            "$25,000",NA))

#recode knowledge
multi2_reshape$scenknow1 <- ifelse(multi2_reshape$scenknow1=="400 before and ${e://Field/numbertotal2} after",
                                   "1","0")
multi2_reshape$scenknow2 <- ifelse(multi2_reshape$scenknow2=="${e://Field/price}%",
                                   "1","0")

multi2_reshape$scenknow_all <- ifelse(multi2_reshape$scenknow1=="1" & multi2_reshape$scenknow2=="1",
                                      "High Knowledge",
                                      ifelse((multi2_reshape$scenknow1=="1" & multi2_reshape$scenknow2=="0") |
                                               (multi2_reshape$scenknow2=="1" & multi2_reshape$scenknow1=="0"), "Medium Knowledge",
                                             "Low Knowledge"))

multi2_reshape$priceafter <- as.factor(as.character(multi2_reshape$priceafter))
multi2_reshape$priceafter <- relevel(multi2_reshape$priceafter, ref = c("Same price"))

#recode into factors
multi2_reshape$numhsafter <- as.factor(as.character(multi2_reshape$numhsafter))
multi2_reshape$numlsafter <- as.factor(as.character(multi2_reshape$numlsafter))
multi2_reshape$product <- as.factor(as.character(multi2_reshape$product))
multi2_reshape$wagehsafter <- as.factor(as.character(multi2_reshape$wagehsafter))
multi2_reshape$wagelsafter <- as.factor(as.character(multi2_reshape$wagelsafter))
multi2_reshape$screen1 <- as.factor(as.character(multi2_reshape$screen1))
multi2_reshape$scenknow_all <- as.factor(as.character(multi2_reshape$scenknow_all))
multi2_reshape$scenknow_all <- relevel(multi2_reshape$scenknow_all, ref = c("Low Knowledge"))
multi2_reshape$self_know_automation_cat <- relevel(multi2_reshape$self_know_automation_cat, ref = c("Low Self-Know"))
#labels
attr(multi2_reshape$numhsafter, "label") <- "Number of High Skilled - After"
attr(multi2_reshape$numlsafter, "label") <- "Number of Low Skilled - After"
attr(multi2_reshape$wagelsafter, "label") <- "Wage of Low Skilled - After"
attr(multi2_reshape$wagehsafter, "label") <- "Wage of High Skilled - After"
attr(multi2_reshape$priceafter, "label") <- "Price - After"
attr(multi2_reshape$product, "label") <- "Product"

#### models #####

#marginal means for DV fairness

pdf(file = "figure3.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69)# The height of the plot in inches
# figure 3
plot(cregg::mm(multi2_reshape, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                 priceafter + product, id= ~ response_id), size=3) + theme(text = element_text(size = 15))  
dev.off()


mm <- cregg::mm(multi2_reshape, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                  priceafter + product, id= ~ response_id) 

mm$outcome <- "fairness"

# amce for fairness

# figure 4

pdf(file = "figure4.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69)# The height of the plot in inches
plot(cj(multi2_reshape, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id),
     size=3) + theme(text = element_text(size = 15))  
dev.off()
cj <- cj(multi2_reshape, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
           priceafter + product, id= ~ response_id)

cj$outcome <- "fairness"

#### Subgroup analysis ######
## fairness ###
## by objective knowledge

#marginal means
mm_by <- cregg::cj(multi2_reshape, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                     priceafter + product + college, id= ~ response_id, estimate = "mm", by = ~scenknow_all)
plot(mm_by, group = "scenknow_all", vline = 0.0, size=3) + theme(text = element_text(size = 15)) 

#amces
amces <- cregg::cj(multi2_reshape, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                     priceafter + product + college, id= ~ response_id, estimate = "amce", 
                   by = ~scenknow_all)
#differences in amces
diff_amces <- cregg::cj(multi2_reshape, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                          priceafter + product + college, id= ~ response_id,estimate = "amce_diff", 
                        by = ~scenknow_all)
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
# only show high and low knowledge for main text figure 5
diff_amces <- filter(diff_amces, BY=="High Knowledge - Low Knowledge")
amces <- filter(amces, BY!="Medium Knowledge")
amces$BY <- droplevels(amces$BY)
levels(amces$BY) <- c("High Know", "Low Know")
diff_amces$BY[diff_amces$BY == "High Knowledge - Low Knowledge"] <- "High - Low" 

pdf(file = "figure5.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces,diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 18))
dev.off()
## subjective self-knowledge
# fairness

#amces 
amces <- cj(multi2_reshape, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product, id= ~ response_id, estimate = "amce", 
            by = ~self_know_automation_cat)
#differences in amces
diff_amces <- cj(multi2_reshape, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product, id= ~ response_id,estimate = "amce_diff", 
                 by = ~self_know_automation_cat)
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 

# high and low self-knowledge only for main text plot figure 6
diff_amces <- filter(diff_amces, BY=="High Self-Know - Low Self-Know")
amces <- filter(amces, BY!="Medium Self-Know")
amces$BY <- droplevels(amces$BY)
levels(amces$BY) <- c("High Self-Know", "Low Self-Know")
diff_amces$BY[diff_amces$BY == "High Self-Know - Low Self-Know"] <- "High - Low" 

pdf(file = "figure6.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces,diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 18)) 
dev.off()

